Ergodic directional switching in mobile insect groups 
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Abstract 

We obtain a Fokker-Planck equation describing experimental data on the collective motion of 
locusts. The noise is of internal origin and due to the discrete character and finite number of 
constituents of the swarm. The stationary probability distribution shows a rich phenomenology 
including non-monotonic behavior of several order/disorder transition indicators in noise intensity. 
This complex behavior arises naturally as a result of the randomness in the system. Its counterin- 
tuitive character challenges standard interpretations of noise induced transitions and calls for an 
extension of this theory in order to capture the behavior of certain classes of biologically motivated 
models. Our results suggest that the collective switches of the group's direction of motion might 
be due to a random ergodic effect and, as such, they are inherent to group formation. 

PACS numbers: 87.23.Cc, 05.40.-a, 05.65.+b, 87.10.Mn 
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I. INTRODUCTION 



Emergence can be defined as the appearance of rich structures on a large scale resulting 
from a multiplicity of simple interactions at a considerably smaller scale. Collective ani- 
mal motion is a paradigmatic example of such an emergent phenomenon. Depending on 
the species, there may exist different hierarchical levels that determine how collective dis- 
placements are realized. For example, in primate groups, an individual's dominance status 
can affect its role in initializing collective movement. In the case of swarming locusts no 
such hierarchies are present; the ability of each individual to guide the band appears to be 
distributed relatively evenly throughout the insect group. Herein we will concentrate on 
groups of wingless locust nymphs which form marching bands rather than flying swarms [1] . 
The onset of collective motion in locusts was experimentally demonstrated in [I], where it 
was shown that sufficiently large insect densities placed in a ring-shaped arena gave rise to 
a coherent displacement of the band. Low densities were characterized by random dispersal 
of the individuals, while for intermediate densities the coherent motion was interrupted by 
sudden changes of direction (hereafter referred to as "switches"). This phenomenology was 
partially rationalized by means of an adapted model based on that of Czirok et al. [2J, who 
formulate a paradigmatic model for collective animal behavior in one dimension. In their 
original model the position, Xi, and velocity, Uj, of locust i are evolved using the following 
two rules, identical for each individual, i = 1, . . . , N, 

Xi(t + 1) = Xi{t) + V Ui(t), 

Ui(t + 1) = G(ui)+&, 

where N is the total number of locusts. Here Ui is the mean of the nondimensionalised 
velocities of locusts within a certain radius, R, of the position, X{, of locust i. The function 
G is such that G{u) = (1 + K)^ 1 ^ + K sgn(u)} for a positive constant K, where sgn(w) 
denotes the sign of u. The role of G is to adjust the average nondimensionalised velocity 
perceived by each particle towards unity. Vq is a constant associated with the chosen time 
scale and is a random number drawn from the uniform distribution in [—77/2,77/2]. The 
adapted version of the model used in pQ to model the movements of locust nymphs in a 
quasi-one-dimensional arena takes the form 
dx ■ 

= m, dui = [G(ui) - Ui\dt + ftdWi, for i = 1, . . . , N, (1) 
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where dWi denotes the increments of independent Wiener processes, (3\ is a positive constant 
describing the amplitude of the noise and the function G is as above. 

A biologically motivated refinement of the model described by Eq. ([I]) was given in [3J, 
where it was postulated that individual locusts increase the randomness of their movements 
in response to a loss of group alignment. This behavior is the result of a particular mul- 
tiplicative form of the noise term (see Eq. (J2|), as opposed to the additive noise in Eq. 
([I]); this characteristic was shown to increase the coherence of the group motion and to 
reduce the frequency of direction switches [3J • The key point in the analysis performed in [3J 
was the estimation of coefficients of an effective Fokker-Planck equation (FPE) [4] , which is 
written in terms of a macroscopic (low-dimensional) observable [5], the average velocity of 
the marching group, derived directly from the experimental data. In the present work we 
approximate the drift and diffusion coefficients of the effective FPE by analytical functions. 
This permits a more thorough analysis and fosters further understanding of collective dy- 
namics of locusts. In addition we compare our results with those of Eq. 0, and discuss the 
disparities between the two models. 

II. THE MODEL 

Coarse-grained analysis [4] allows us to obtain an effective FPE describing the collective 
behavior of the locusts at the macroscopic level. By using this coarse-graining technique (see 
[3j) we were able to extract the coefficients of the assumed underlying FPE describing the 
alignment of the locusts from the experimental data presented in [1] . This approach enables 
us to reduce our system - comprising a large number of degrees of freedom - to a single 
collective variable, u, (referred to variably, hereafter, as 'alignment' or 'average velocity') 
which characterizes the system's macroscopic behavior. The proposed FPE has a simple 
form and it can be expressed as 

+ §cU(i-u 2 )n (2) 

for the probability P(u, t)dudt of finding the system with an average velocity in the interval 
(u, u + du) during the time interval (t,t + d£); note that the experimental situation in [1] 
is quasi-one- dimensional, allowing the use of a one-dimensional FPE [3j. We note that this 



d t P = -a 2 d u 



u 



u — 
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FPE corresponds to the following Langevin equation for the average velocity 



du = a 2 (u--— 2 )dt+J^VT^dW, 




(3) 



where dW denotes the increments of a Wiener process and the multiplicative noise is in- 
terpreted in the Ito sense, as prescribed by the experimentally obtained FPE (|2]). In these 
equations the average velocity u is dimensionless and takes its values in the interval [—1, 1]. 
The values such that \u\ = 1 characterize the ideal situation in which all locusts march in 
perfect coherence; the sign determines the direction. Of course, formally substituting \u\ = 1 
in Eq. ^ produces a divergence in the drift, so for practical reasons one has to assume that 
coherent motion implies \u\ < 1 rather than a strict equality. The value u = characterizes 
a total disorder; realistic values of the average velocity lie between these two extreme cases. 
Although the average velocity is a dimensionless quantity in Eqs. ^ and ([3]), time is not. 
Consequently a 2 and f3 2 have the dimensions of time -1 . We shall estimate in the following 
their numerical values using the experimental data from pQ and express them in units of 
seconds -1 . 

The proposed FPE ^ describing the alignment is a reasonably accurate approximation 
to the unknown FPE assumed to underly the motion of the locusts, which captures their 
experimental swarming behavior. It should be noted that such an equation can only be 
obtained if the system being studied is amenable to this sort of reduction. 

For asymptotically large values of N, parameter 1 (in equation (|2])) denotes the order of 
magnitude of the relaxation time characterizing how long it takes the entire group to become 
ordered when starting from a disordered configuration, and N fi^ 1 indicates the order of 
magnitude of the characteristic time over which the fluctuations of the mean velocity develop. 
For the range of experimentally considered locust numbers (5 < N < 40) the observed values 
of a 2 and f3 2 are approximately constant while we expect the presence of a boundary layer 
for smaller values of N. Since our results in [3] are rather noisy our goal is to fit the order 
of magnitude of the model parameters instead of attempting to obtain precise estimates. 
Comparing the proposed analytical coefficients of Eq. ^ to those obtained in [3] , from the 
experimental data in PQ, we obtain f3 2 /a 2 = 2.4 ± 1.7. Employing the mean switching time 
measurements in [3] we find a 2 = (6.65 ± 2.63)10 -4 s -1 and (3 2 = (1.62 ± 0.52)10 -3 s -1 . 

The FPE corresponding to Eq. (TO) can be obtained as a mean-field approximation, 




(4) 
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where a\ = K/(l + K) and K is defined as for the function G in Eq. ([I]). The stationary 
solution of the FPE Q can be derived as follows: 



2tt^i I 2/3 
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The values of the two maxima of this stationary probability distribution (SPD) M max = ±1 
and the minimum u m ; n = are independent of the parameter values. This type of system 
has been considered many times in the literature [6j [7] , and we include it here simply for 
completeness and comparison with the refined model (|3]). 

In the absence of sources and sinks of probability, we can also derive the SPD of the 
experimentally motivated FPE (J2|: 

Na 2 /(2f3 2 



P S ( U ) =M{\ _ M 2)-l-AW/3 2 
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where M~ l = /^(l — u 2 )~ l ~ Na2 ^ 2 exp [— A^a2/[2/?2(l — w 2 )]] du is the inverse of the normal- 
ization constant. This SPD is bounded, compactly supported in [—1, 1] and bimodal for all 
values of the parameters. 

Noise induced transitions have been studied traditionally by means of the dynamics of the 
extrema of the SPD [8] . For the biologically motivated FPE ^ the SPD shows one minimum 
always located at w min = 0, and two maxima at u max = ±a/«2 + Zfo/N / \/2a 2 + 2/3 2 /iV, 
which exist for all parameter values. One immediately notes |w ma x| £ (l/\/2, 1), a fact 
related to the shape of the "deterministic potential" (the potential in the absence of noise), 
which is the negative integral of the drift coefficient, 



o 



s 3 



V(u) = - / a 2 [s - ds = -a 2 



u 2 + hn(l-u 2 ) 



(7) 



This potential is bistable with one maximum located at the origin and two minima at 
±1 / a/2 independent of the parameter values. For increasing noise intensity the probability 
maxima of the SPD ^ (corresponding to the biologically motivated FPE ([2])) separate 
from the deterministic potential minima ±l/v2 and approach the boundary points ±1. 
These facets of the SPD, when considered in the context of the classical theory of noise 
induced transitions, imply that the system is becoming ordered j8]: the SPD maxima, 
representing the states in which the system will most likely be found, are further apart 
and thus there is a clearer differentiation among those states. However, the experimental 



evidence, based on switching times which decrease as the noise magnitude increases, reveals 
that the system becomes disordered [Tj. This indicates that for complex systems, restricting 
the characterization of the dynamics to observations of the evolution of the extrema may 
not be adequate in some experimentally motivated situations. Herein we will try to carry 
out a more complete characterization. 



III. BARRIER HEIGHT 



Another indicator of order /disorder is the barrier height of the effective potential. For 
the model given by Eq. (1) the barrier height decreases monotonically as the noise intensity 
increases as can be seen from Eq. (|8l 



a>iN ( , , 1 2 \ a,, _ / \ / \ a i N 



Vcs(u) = — ^— y\u\ - -u j , AV cff = V cff (u min ) - V cff (w max ) = (8) 

where V e s is the effective potential and AV e ff the corresponding barrier height. 
The effective potential for our revised model (Eq. (l2l) is given as 



VMu) ^ - ln[P,H = + (l + g ln(l - A (9) 

and the corresponding barrier height is 

AV cff = V cff ( Mmin ) - V cff ( Mmax ) = -1 - ^ + (l + ^\ In (2 + . (10) 



min ~ 



As a function of noise intensity the barrier height exhibits a minimum at [/^/(-/Va^)] 
0.76. This means that for (32/{Na 2 ) < \fi%l (Na 2 )] m i n (sub-threshold) the barrier height di- 
minishes for stronger noise, but for /3 2 /(Aa; 2 ) > [/^/(iVa^jmin (super-threshold) it increases 
as the noise strength grows. Indeed, AV e ff ~ [ln(2) — 1/2] (Na 2 / ^2) when Na 2 / (3 2 — > 00 and 
AV c ff ~ — \ri(Na 2 / (3 2 ) when Na 2 /(3 2 — > 0. This suggests that, while increased noise causes 
the system to become more disordered for sub-threshold noise intensities, super-threshold 
intensities might cause the system to become more ordered as the noise grows. In short, 
the 'barrier height' order parameter shows a clear non-monotonicity when considered as a 
function of noise strength. This appears like a counterintuitive reentrant behavior, where 
the noise can have an ordering effect for supercritical intensities [8] . Although this behavior 
is interesting in itself, it is not biologically relevant, as it requires a number of individuals 
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FIG. 1: (Color online) Panels (a) and |(b)| show the profile of the SPD, P s {u), from Eq. ^ plot- 
ted against the normalized mean velocity u for a varying noise intensity f5i/{Nai2). In panel [(aj 
the solid red line represents f32/(Na2) = 0.05, the dashed green line represents ^/(-^a^) = 0.15 
and the dotted blue line represents f32/(Nct2) = 1.5. In panel [(b)] the solid red line represents 
/^/(-/Va^) = 1.5, the dashed green line represents $2/^0,2) = 3.0 and the dotted blue line repre- 



sents /^/(iVa^) = 5.5. Panel (c) displays the second moment of the revised model, centered at the 
origin, So, versus noise strength /^/(-^Q^) and panel [(d)] shows the second moment at a maximum, 
S m , versus noise strength ^2/ '(iVa^)- The minima are attained for /^/(-^Va^) ~ 0.27 (So) and for 
f3 2 /(Na 2 ) « 0.12 (S m ). 

N w 3, beyond the validity of the model. Both characteristics of the SPD ^ of our re- 
fined model, displacement of the location of the maxima and non-monotonic variation of the 



barrier height, can be seen in Fig. 1 (a) (b) 



Let us note that noise induced non-equilibrium phase transitions [HI [H)] as well as other 
noise mediated ordering phenomena [11] have been exhaustively studied in the literature. 
This includes the detailed study of reentrant transitions [12J. However, most of these 
approaches have assumed multiplicative noise interpreted in the Stratonovich sense. The 
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Stratonovich interpretation is able to destabilize a deterministically stable state and in 
this way produce phase transitions or other noise induced phenomena. The key technical 
point in these cases is the appearance of a systematic contribution to the deterministic 
dynamics coming from the noise term, the so-called "Stratonovich drift" [8]. Therefore noise 
interpretation plays a fundamental role in the development of these types of phenomena: 
in particular, many noise induced phenomena are not possible if the noise interpretation 
is that of Ito. In this respect, our results are fundamentally different as we only consider 
the Ito interpretation for our Langevin equation with multiplicative noise ([3]). We also note 
that noise induced phase transitions which are independent of the noise interpretation have 
also been studied [13] . but in much less detail. The mechanisms leading to these phase 
transitions are based on the bifurcation of the minima of an effective potential due to noise 
and, as such, constitute the natural extension of noise induced transitions [8] to spatially 
extended systems. Our results relate to a zero dimensional system as do those in [8], but 
they are significantly different as the model defined by Eq. ^ does not describe this kind 
of bifurcation. 



IV. MEAN SWITCHING TIME 



We can further explore the model properties by considering the mean switching time, 
T(u), defined as the first time, on average, that the alignment of the system, u, initialized 
such that — 1 < u < 0, reaches the origin (u — 0). For our revised model the moments of 
the switching time distribution are given, recursively, by the solution of the equation 

u 3 \ „ „ (3 2 



CX 2 U 



U 



u )d uu T n = -n% 



n— 1 1 



(11) 



subject to the boundary conditions T(0) = and T'(— 1) = 0, where T n is the n th moment, 
correspondingly, T = T\ is the mean switching time and Tq = 1. The second boundary 
condition represents zero probability flux through u — — 1. Eq. (11 ) is directly derived from 



the FPE using methods from [TJj. The solution to this equation for n = 1 is 



T(u) 
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This expression appears complicated, but one can derive its asymptotic expansion for large 
values of Na 2 /(3 2 (which implies large N as a 2 and f3 2 are approximately constant) by means 



i— f \ N012/ @2 

of a Kramers-like approximation |15j . It has the simple form T [— l/\/2) ~ ( ^ J > 
which reveals a pure exponential growth in the inverse noise intensity N 0L2/ P2 for asymp- 
totically large values. Further moments of the switching time distribution can be calculated 



from Eq. (11) for n > 1. In the limit N — > 00 one finds the relation T n = n\T n . This 
relationship implies, in turn, that the switching process is a Poisson process. 

We can also compute the first passage time for the model Eq. Q. In this case we solve the 
equation a\[sgn(u) —u}d u T + ^d uu T = —1, subject to the boundary conditions T(0) = and 
T'{— 00) = 0, where the latter condition is the analogue of the previous zero flux condition 
adapted for an SPD with infinite support. We find 







r(tl) = ,/ 2^' 6XP 



2 — erfc 



' OL\N , 



dv. 



which also behaves exponentially in N for large values of iV but this time with an N de- 
pendent prefactor (see Supplementary Information of [3]). The relation between these two 
mean switching times (the model Eq. (II]) and that of the revised model [3J) is extensively 
discussed in [3], so we will not reproduce the discussion here. 

Now we compare the theoretical results with the exponential fitting we have performed 
on experimental data from [1] for both the first and second moments of the switching time 
distribution. The data are insufficient for us to reliably obtain any moments higher than 
the second. For the mean switching time T and second moment T2 we found 



T = (970 ± 120) exp[(0.045 ± 0.007)iV] s, ^/%/2 = (1300 ± 190) exp[(0.041 ± 0.008)iV] s. 

According to the relation T n = n\T n for the moments of the exponential distribution, these 
two values should be the same if the switching process were Poissonian. Note that the ex- 
ponential growth is the same for both (within errors), while the prefactor is larger for the 
second moment. This suggests that the switching process is Poissonian for large N, that 
is, the probability distribution for the switching events is V — T -1 exp(— t/T). For small 
values of N the behavior is more stochastic, as signaled by the larger prefactor of the second 
moment (when N is small the prefactor dominates over the exponential). If the switching 
process is Poissonian then this has a series of consequences concerning predictability: the 
standard deviation being equal to the mean implies a 100% error in predictions. Further- 
more, switching events are uncorrelated and the distribution tail falls off exponentially for 
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long times. This allows for a higher probability of rare events than would be allowed by a 
Gaussian tail. This also implies that the switching process is Markovian, as predicted by the 
FPE. This can be seen from the double-welled FPE ^ in the large N limit. After a short 
time the system relaxes to one potential minimum where it stays an exponentially long time 
until the switch occurs. Since practically all switches start at the minimum this erases the 
memory and the Markov property is recovered. The verification of this theoretical prediction 
by the experimental data suggests that no important correlations have been suppressed in 
the coarse-grained computation in [3], and that this method, and the FPE Q, are suitable 
to describe the locust dynamics exhibited by the experimental data. 



V. SECOND MOMENTS 



Another indicator of the stochastic properties of the system is the second moment, which 
measures the spread of the mean velocity, u, with respect to some reference value. We 
consider two variants, one centered at the origin So = f_, u 2 P s (u)du, and one centered at 
one of the maxima of the probability distribution S m = J_ t (u — u max ) 2 P s (u)du. Of course, 
the value of S m is the same for both maxima as a consequence of the symmetry of the system. 
These integrals have been computed numerically and are represented in Fig. l|[c)| centered 
at the origin, and l|[d) centered at a maximum. Both show non-monotonic behavior in noise 
intensity, but attain their minima for different values of the noise amplitude. This non- 
monotonic behavior, as well as the behavior of the effective barrier height, are not reflected in 
the relationship between mean switching time and the size of the noise parameter (N0L2/ fa)'- 
the mean switching time grows monotonically with noise amplitude. For comparison we note 
that both moments So and S m grow monotonically with the inverse noise intensity in the 
model given by Eq. ([!]); in this case they are 



>->0 — 1 H H 7 — / n — , >->m — * H 77 H 7 — ; x — • U^J 

aiN ! + - f (\/W) ^ 1 + erf (vW 

There is another feature of the second moments of the revised model |3j, in addition to 
the non-monotonic behavior, that reveals new characteristics of the collective motion of 
locusts not reflected by the model Eq. (II]). In this model a reduction in the number of 
individuals increases the values of both second moments. In the stronger noise situation 
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the probability distribution tails grow, which implies that there are more individuals with 
a higher (absolute value) velocity. In our case the probability is compactly supported in 
[—1,1], as a consequence of the biological fact that the propagation cannot be better than 
perfect. For realistic values of the parameters the system is in the weak noise regime (see Fig. 



1 c) ). This means that the second moment centered at the origin decreases for a decreasing 



number of locusts, exactly the opposite trend to that of the model Eq. ([I]). The reason is that 
the probability of finding the system in the neighborhood of u = grows considerably for 
stronger noise (as reflected by the decreasing barrier height), largely compensating for the 
drift of the maxima towards the boundaries of the support of the SPD. The experimentally 
derived value of P 2 /Na2 = 0.12 ± 0.08 for N = 20 agrees with the minimizing value of the 
second moment centered at a maximum, (3 2 /Na 2 « 0.12 (see Fig. |l|(d)| ). This implies that 
its behavior is not very sensitive to small changes in the number of locusts. 



VI. CONCLUSIONS 



We have seen that the FPE obtained from the coarse-grained analysis of experimental 
data on the movement of locusts shows an interesting phenomenology. Different indicators of 
order /disorder may vary non-monotonically with noise intensity, possibly in a contradictory 
manner. These findings reveal that these indicators might not be suitable for the biologically 
motivated models studied in this paper. We have also shown that the direction switches 
are independently distributed for large numbers of individuals. This makes them almost 
unpredictable from a practical viewpoint. It seems that directional switches are produced by 
an accumulation of errors (made by the locusts when trying to adapt their velocity to that of 
their neighbors) that ordinarily interfere and cancel each other out but, over exponentially 
long times, have the possibility of accumulating and producing a switch. According to 
the results presented here, specifically the confirmation of the Poissonian character of the 
switching events, it seems possible that directional switches are produced as a consequence 
of the ergodic random evolution of the system. We note the similarity of this process with 
Ising model ergodic magnetization changes [TU]. Indeed, the model of Eq. ([I]) can be thought 
of as an Ising model with moving spins. It seems that the ergodic nature of the finite size 
Ising model is preserved despite introducing movement of the spins. More importantly it 
seems that this is a plausible explanation, in the absence of external stimuli, of the sudden 
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changes of direction observed in animal groups. 
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